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1. Motivation and Objectives 

Transition to turbulence in aerospace applications usually occurs in a strongly dis- 
turbed environment. For instance, the effects of free-stream turbulence, roughness 
and obstacles in the boundary layer strongly influence transition. Proper unde 
standing of the mechanisms leading to transition is crucial in the design of aircraft 
wings and gas turbine blades, because lift, drag and heat transfer strongly depend 
on the state of the boundary layer, laminar or turbulent. Unfortunately, most of 
the transition research, both theoretical and experimental, has focused on natural 
transition. Many practical flows, however, defy any theoretical analysis and are ex- 
tremely difficult to measure. Morkovin 5 introduced in Ms review paper the concept 
of bypass transition as those forms of transition which bypass the known mecha- 
nisms of Unear and non-Unear transition theories and are currently not understood 

by experiments. . , . . . . „„ 

In an effort to better understand the mechanisms leading to transition in an 

disturbed environment, experiments are conducted studying sunpler cases, viz the 
effects of free-stream turbulence on transition on a flat plate, Sohn and Resho o 
and Wang et al. 19 . It turns out that these experiments are very difficult to conduct, 
because the generation of free-stream turbulence with sufficiently high fluctuation 
levels and reasonable homogeneity is non trivial. For a discussion see Morkovin . 
Serious problems also appear due to the fact that at Mgh Reynolds numbers the 
boundary layers are very thin, especially in the nose region of the plate where the 
transition occurs, which makes the use of very small probes necessary. 

The effects of free-stream turbulence on transition are the subject of this re- 
search and are especially important in a gas turbine environment where turbu- 
lence intensities are measured between 5 and 20%, Wang et al. 19 . Due to the fact 
that the Reynolds number for turbine blades is considerably lower than for aircraft 
wings, generally a larger portion of the blade will be in a laminar-transitional state. 
Turner 15 shows that the effect of free-stream turbulence on transition significantly 
increases when the free-stream turbulence levels become larger than 5% arid is ac- 
companied with a significant increase in heat transfer. Recently Rai and Mom 
presented a direct numerical simulation of transition to turbulence on a flat plate 
in a free-stream with Mach number .1 and turbulence levels at the leading edge 
of about 2 75%. Direct numerical simulations offer a unique opportunity to study 
specific phenomena, while excluding disturbances from other sources The com- 
putations from Rai and Moin show some impressive results, especially regarding 
intermittency and turbulent spots. Their numerical simulation, however, has the 
same problem as with most of the experiments, namely a very low Mach number, 
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while many applications operate in the transonic regime. Due the nature of their 
numerical scheme, a non-conservation formulation of the Navier-Stokes equations, 
it is a non-trivial extension to compute flow fields in the transonic regime. 

This project aims at better understanding the effects of large free-stream tur- 
bulence in compressible boundary layers at Mach numbers both in the subsonic 
and transonic regime using direct numerical simulations. The present project aims 
at computing the flow over a flat plate and curved surface. This research will 
provide data which can be used to clarify mechanism? leading to transition in an 
environment with high free stream turbulence. This information is useful for the 
development of turbulence models, which are of great importance for CFD appli- 
cations, and are currently unreliable for more complex flows, such as transitional 
flows. 

2. Accomplishments 

Direct simulations of transition in compressible flows with both shocks and bound- 
ary layers requires an extremely accurate and efficient scheme. Several conflicting 
requirements present a serious challenge which cannot be met by existing numerical 
schemes: 

• The small grid spacing in the boundary layer makes an implicit scheme necessary, 
because an explicit scheme would have a severe time step limitation. Implicit 
schemes usually are not time accurate and rather dissipative. • 

• Higher order accurate schemes are necessary but higher order accurate schemes 
generally do not give non-oscillatory solutions around discontinuities, such as 
shocks. Many of the popular non-oscillatory shock capturing schemes, such as 
TVD (Total Variation Diminishing) methods, are only first order accurate in 
multi-dimensional flows and even in one-dimension they reduce to first order at 
non-sonic local extrema. 

In order to satisfy these conflicting requirements a significant effort has been made 
to improve and combine several successful numerical schemes. A fully implicit and 
time accurate code for the solution of the three-dimensional compressible Navier- 
Stokes equations in general geometries has been written and tested. Higher order 
accuracy and shock capturing are implemented using an Essentially Non-Oscillatory 
(ENO) scheme. Time accuracy is obtained using a Newton method. 

In the next section a brief description of the numerical scheme will be given 
followed by the discussion of a series of tests aimed at validating the code. 

2.1 Numerical Scheme 

The compressible Navier-Stokes equations are solved using a finite volume method ! 
A detailed discussion of finite volume and difference methods can be found in 
Vinokur 18 . The integral formulation of the Navier-Stokes equations, assuming all 
variables are continuous in time, is given by: 

— f U dV+<f n • TdS = 0 (2.1.1) 

ot JV(t) Js(t) 
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Here V(t) and surface S(() are the volume and outer surface of the domain n 
and n an outward unit normal vector at 5. The vector U represents the conserved 
variables: (p, pu, pv, pro, ef , with p the density u and " 

and e the total energy. The tensor 7 is defined as T= £+ V, with £ the invisc 
contribution defined as: 
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The shear stress tensor % with components (ti,t 2 ,t 3 ) is defined as: 
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and the heat flux q as: 


kVT 


(2.1.5) 

q ~ (7 - l)ReMlPr 

The variables p, T, //, A and /c represent the pressure, temperature, fi^taiid 
second viscosity coefficient and thermal conductivity, respectively. The coefficients 
Re, Moo, and Pr are the Reynolds, Mach, and Prandtl numbers. A1 variables are 
non-dimensionalized using free-stream variables and a characteristic length L. 

The Navier-Stokes equations are solved using a finite volume method because we 
seek a weak solution in order to capture shocks in high Reynolds number flows. 
The finite volume method is also the most natural way to satisfy the conservation 
properties of the differential equations. After subdividing the vo ume V into a set of 
disjunct cells we obtain the finite volume discretization for a cell with index i,], • 




d_ 

dt x 3 ,J '~ -w * 2'-- - *• - (2 16) 

where a barred quantity with index i,j,k is an average of the unbarred quantity 
over the cell with index i,j,k and indices i±±, j and fc ± * refer to values a 
the cell faces. The fluxes F* at the cell faces are defined as: 

F { = S l T (2.1.7) 

with S l the cell face in the direction i. The computation of the cell face S 1 and 
volume V has to be done with great care in order satisfy the geometric conservation 

law, for details see Vinokur 19 : 
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Flux Approximation 

The crucial part in the development of a finite volume method is the approxi- 
mation of the fluxes at the cell faces. The flux F 1 j is computed using the Osher 

l-h 2 

approximate Riemann solver. The first order accurate conservative flux is given by: 




+ F], 





( 2 . 1 . 8 ) 


with equivalent expressions for the other two directions. The integral is computed 
along a path in phase space, connecting the points with index ( i,j,k ) and (t+1, j, k). 
Along each subpath a Riemann problem is solved, which is used to determine the 
intermediate states. In this way exact expressions can be derived for the path inte- 
grals. More details about the implementation of the Osher scheme can be found in 
Osher and Solomon 6 , Osher and Chakravarthy 7 , Chakravarthy and Osher 1 and Rai 
and Chakravarthy 10 . The Osher approximate Riemann solver is the most accurate 
approximate Riemann solver and satisfies the entropy condition, contrary to the 
Roe approximate Riemann solver which needs an entropy fix to eliminate steady 
expansion shocks. The Osher scheme captures steady shocks in at most two points. 

The most important reason for the choice of the Osher scheme, however, has 
been its low numerical dissipation in boundary layers, Koren 4 . Most schemes for 
the Euler equations are very dissipative in the boundary layer and not well suited 
for direct numerical simulations. In earlier work, Van der Vegt 17 , modifications to 
flux vector splitting schemes were discussed to alleviate this problem, but although 
significant improvement was achieved on steady laminar boundary layers, it was 
not possible to reach accuracy levels necessary for direct simulations. 

Higher order spatial accuracy 

Direct simulations require a high accuracy which cannot be achieved with stan- 
dard second order schemes. It is fairly straightforward to derive higher order accu- 
rate finite difference schemes, but shock capturing then will not be possible. The 
development of higher order accurate, multi-dimensional finite volume schemes, 
capable of shock capturing still is an area of active research, and has been an im- 
portant subject in this project. A significant effort has been made to combine the 
Osher approximate Riemann solver, discussed in the previous section, with an ENO 
scheme. Results of this work are described in Van der Vegt 17 , where the different 
ENO methods are discussed and results of various tests are discussed. 

Higher order accuracy of a finite volume method can be defined is various ways. 
One approach is to define higher order accuracy with respect to the cell averaged 
values. This resembles most closely the finite volume description, which gives equa- 
tions for the cell averaged values. Another definition of higher order accuracy uses 
the point values at the cell center, which is used in conservative higher order finite 
difference methods. Both approaches are being used. For subsonic flows currently 
the fifth order scheme, developed by Rai 9 , is used, which is based on a Taylor series 
expansion of the flux vector along the lines presented by Osher and Chakravarthy 8 . 
This method is a conservative finite difference scheme. It has the benefit that it is 
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simple to implement in more dimensions, but does not allow shock capturing. Os- 
her and Chakravarthy 8 demonstrated how to make these schemes TVD and allow 
shock capturing, but they are not very useful and only first order accurate global y. 
The scheme therefore is used in its unlimited form, limiting its application to flows 
without discontinuities such as shocks. The scheme also is rather expensive and 

work is in progress to improve its efficiency. 

For flows with shocks research has been carried out to develop higher order ac- 
curate ENO schemes. ENO schemes use an adaptive stencil where a searching 
algorithm tries to find that part of the flow field surrounding a cell which is the 
smoothest. Then a conservative, higher order accurate interpolation method is used 
to ’’reconstruct” the point values from the cell averaged values. Due to the fact that 
the interpolation process only uses data from the smooth part of the flow field nu- 
merical oscillations will be minimized. In this way uniform higher order accuracy 
can be obtained. The first ENO methods were developed by Harten etal and 
later modifications were proposed by Shu and Osher • . In Van der Vegt t e 
different methods were compared and it was found that the ENO scheme, using 
primitive function reconstruction from cell averaged variables with the Cauc y- 
Kowalewski procedure for time integration combined with the Osher approximate 
Riemann solver, was the most accurate and robust. In one dimension it has been 
successfully used up to fifth order accuracy, but due to the fact that in multi- 
dimensional flows currently dimension splitting is used, its accuracy is limited to 
second order in more than one dimension. Work to extend this scheme to higher 
order accuracy in multi-dimensional flows is in progress. 


Time integration 

Due to the very small gridspacing necessary at the wall and in critical layers 
explicit time integration would result in a serious time step limitation. To alleviate 
this problem implicit time integration has to be used, but most implicit time inte- 
gration schemes make assumptions in the implicit part which reduce or eliminate 
time accuracy. The development of implicit and time accurate numerical schemes 
therefore has been a significant part of this research. Time accuracy is obtaine 
using the Newton method discussed in Rai 9 , which solves the non-linear system of 
equations in the implicit time integration scheme using a Newton method. Rai uses 
this method also to reduce the error caused by approximate factorization. We do 
not use approximate factorization but solve the whole matrix system iteratively, see 
Van der Vegt 16 , and need the Newton scheme only to reduce the error due to the 
time linearization. This iterative scheme also has the benefit that it is not neces- 
sary to use an exact linearization of the flux vectors, which can be very difficult 
and time consuming to obtain. First order Steger- Warming flux vector splitting is 
used in the implicit scheme, while a higher order accurate spatial discretization is 
used for the explicit part. At each time step the Newton iteration is performed 
such that the accuracy of the higher order explicit part is maintained. The use of 
an approximate linearization, however, limits the maximum time step and work is 
in progress to evaluate if a more accurate linearization would improve the perfor- 
mance and robustness of the scheme. Especially at high Mach numbers there still 
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are convergence problems for large Courant numbers when the scheme is used to 
obtain steady state solutions. 

2.2 Results and Discussion 

Several computations have been carried out to test the code and the various 
numerical schemes. In order to test the Osher approximate Riemann solver and 
the ENO schemes a large number of shock tube calculations have been carried out, 
a detailed description of this work can be found in Van der Vegt 17 . The different 
ENO schemes tested were the ENO method from Harten et al. 3 , using primitive 
function reconstruction and either Runge-Kutta time integration or the Cauchy- 
Kowalewski procedure, and the Shu and Osher flux-ENO scheme. In all cases the 
Osher approximate Riemann solver was used and the effect of the ordering of the 
eigenvalues, viz. natural and reversed ordering, has been investigated. Four cases 
with different difficulties were tested, see Table 1. The performance of the different 
schemes was reasonable in most cases, but it turned out that the ENO scheme 
with primitive function reconstruction and the Cauchy-Kowalewski procedure for 
time integration (ENO-CK) was the most accurate and robust. Some of the results 
obtained with this method, are shown in Figures 1 and 2. Figure 1 shows a left 
moving shock followed by a contact discontinuity and a right moving expansion 
wave. A difficult problem for the ENO schemes in case A is the fact that in the 
initial stages there are not enough grid points available in the region between the 
shock and contact or shock-shock. The ENO scheme searches for the stencil which 
gives the smoothest part of the flow field around a grid point or cell. In these 
cases there exist in both directions a discontinuity and there are not enough points 
available to build a non-oscillatory higher order reconstruction. This problem exist 
for all ENO schemes but the ENO Cauchy Kowalewski scheme is the least sensitive 
for it. The other methods have mild to strong oscillations in these areas. One way 
to solve this problem is to reduce the accuracy locally till enough grid points are 
available to create a higher order reconstruction, but this is a problem which still 
needs further attention. Case B, Figure 3, shows a left moving expansion wave and 
a right moving contact discontinuity and shock. One of the problems in this case is 
the appearance of a sonic point, which gives a small jump of O(Ax) at first order. 
The shock tube tests showed that it is possible to use a higher order scheme for 
flows with discontinuities, but the convergence of these higher order schemes can 
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A 

15000 

98400 

0 

0 

1378 

4390 

B 

988000 

9930 

0 

0 

2438 

2452 

C 

10000 

100000 

0 

0 

2627 

272 

E 

573 

22300 

2200 

0 

199 

546 


Table 1. Initial Conditions Shock Tube Tests. 
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be at most first order around these discontinuities. Tests of all the ENO schemes 
on smooth solutions showed that they all reached the proper level of accuracy. Test 
are currently underway to check if the ENO schemes give higher order accuracy in 
regions outside discontinuities as they are supposed to. This is an important test 
to see if these methods are capable of shock-turbulence interaction simulations. 

In order to test the shock capturing properties of the code the flow field around 
a circular cylinder at Mach 8 has been computed. Although the flow field was two- 
dimensional the three-dimensional code was used to check if the flow field remained 
exactly two-dimensional and the geometric conservation law was satisfied. Figure 
3. presents contours of the pressure at a spanwise station along the cylinder. e 
solution is the same at all stations. The flow field consists of a strong bow shock 
where the Mach number changes from 8 to about 2.8 behind the shock at the 
symmetry line. Apart from the strong shock another aspect of this case is the 
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Figure 3. Pressure field of flow around a cylinder at Mach 8 


fact that the flow field in the stagnation region ahead of the cylinder is subsonic. 
The sonic line is at about 45 degrees with the flow angle and a smooth transition 
is observed from the subsonic to the supersonic region. This case has also been 
computed by Osher and Chakravarthy 7 and the results compare well. To test the 
ability of the code to simulate transitional flows which is a crucial test before bypass 
transition can be simulated currently computations are done on natural transition 
in a flat plate boundary layer. A comparion is made with the results of Fasel et 


122 





Bypass Transition in Compressible Boundary Layers 



Figure 4. Vertical velocity field in a flat plate boundary layer at Mach 0.08 with 
periodic suction and blowing (vertical direction enlarged 20 times). 


a l. 2 , which computed transition in an incompressible boundary layer. In orde 
make the comparison as accurate as possible a very low Mach number > - M 
was chosen This Mach number is approximately the lower limit for the num 
scheme and despite the fact that the computations are fully implicit ^evere time 
step limitation is imposed by the sound waves. To start the simulation first a steady 
laminar boundary layer is computed, which is a non-trivial task because a ver y. h j6* 
accuracy is needed in this computation. The disturbances at the beginning o 
plate have an amplitude of 10" 4 and transition simulations require a very tow ‘ 
merical dissipation. The disturbances are generated by periodic suction and blowi g 
rrstrip somewhat downstream of the inflow boundary. This is done because there 
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are always numerical disturbances due to the fact that the inflow boundary condi- 
tions are not perfect. The disturbances are generated in a region of the boundary 
layer which is linearly stable. When they move downstream first the transients are 
damped and after they move in the unstable region the Tolmien Schlichting waves, 
which are the most unstable two-dimensional waves, are amplified. In order to ac- 
commodate for the fact that the subsonic outflow boundary conditions, which are 
essentially inviscid, are not perfect in a boundary layer a buffer layer is added to the 
plate to damp as much as possible the reflections coming from this boundary layer. 
This is the same procedure as used by Rai and Moin 11 . The number of grid points 
in this computation is 340 x 82. It turned out that the most efficient way to obtain 
the steady boundary layer was to first start with the second order scheme running 
at a Courant number of 120 using the implicit Euler time integration and to switch 
to the fifth order scheme after most of the transients are damped out. The fifth 
order scheme has a very low dissipation and would otherwise take a long time to 
converge. The maximum CFL number for which the fifth order scheme is stable is 
approximately 160. Ater the steady state was obtained periodic suction and blowing 
were added and the fifth order scheme was used with the Newton time integration 
scheme at a CFL of 60. Figure 4. shows a preliminary result of this computation 
and it clearly shows the gradual build up of the boundary layer instability. The 
results are currently analysed to make a comparison with the incompressible results 
of Fasel et al. 2 . 

3. Future Plans 

Further testing of the code will have to be done for more complicated cases. 
Currently a comparison with the results of stability of a flat plate boundary layer 
at low Mach number with the results from Fasel et al. 2 is being completed. If 
this comparison and equivalent tests for high Mach number boundary layers are 
satisfactory a simulation of bypass transition on a subsonic flat plate boundary 
layer will be made, followed by simulations of a boundary layer on a curved plate. 
Also work has to be continued to develop higher order accurate ENO schemes for 
multi-dimensional flows. This is crucial for direct simulations of transonic flows 
with both shocks and turbulence. 
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